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We present and compare different numerical schemes for the integration of the variational equa- 
tions of autonomous Hamiltonian systems whose kinetic energy is quadratic in the generalized 
momenta and whose potential is a function of the generalized positions. We apply these techniques 
to Hamiltonian systems of various degrees of freedom, and investigate their efficiency in accurately 
reproducing well-known properties of chaos indicators like the Lyapunov Characteristic Exponents 
(LCEs) and the Generalized Alignment Indices (GALIs). We find that the best numerical perfor- 
mance is exhibited by the 'tangent map (TM) method', a scheme based on symplectic integration 
techniques which proves to be optimal in speed and accuracy. According to this method, a sym- 
plectic integrator is used to approximate the solution of the Hamilton's equations of motion by the 
repeated action of a symplectic map S, while the corresponding tangent map TS, is used for the 
integration of the variational equations. A simple and systematic technique to construct TS is also 
presented. 

PACS numbers: 45.10.-b, 05.45.-a, 02.60.Cb 



I. INTRODUCTION 

Numerical integration is very often the only available 
tool for investigating the properties of nonlinear dynam- 
ical systems. Different numerical techniques [ll, Q have 
been developed over the years which permit the fast and 
accurate time evolution of orbits in such systems. 

Of particular interest are the so-called 'symplectic inte- 
grators ' which are numerical methods specifically aimed 
at advancing in time the solution of Hamiltonian sys- 
tems with the aid of symplectic maps (see for example 0, 
Chapt. VI], and references therein). Another challeng- 
ing numerical task in conservative Hamiltonian systems 
is to discriminate between order and chaos. This distinc- 
tion is a delicate issue because regular and chaotic or- 
bits are distributed throughout phase space in very com- 
plicated ways. In order to address the problem several 
methods have been developed, which can be divided into 
two major categories: the ones based on the study of 
the evolution of deviation vectors from a given orbit, like 
the computation of the maximal Lyapunov Characteris- 
tic Exponent (mLCE) xi 0: those relying on the 
analysis of the particular orbit itself, like the frequency 
map analysis of Laskar Q . 

Other chaos detection methods, belonging to the same 
category with the evaluation of the mLCE, are the fast 
Lyapunov indicator (FLI) @ and its variants 0, the 
smaller alignment index (SALI) Q and its generaliza- 
tion, the so-called generalized alignment index (GALI) 
S [13 ) and the mean exponential growth of nearby or- 
bits (MEGNO) The computation of these indicators 
require the numerical integration of the so-called varia- 
tional equations, which govern the time evolution of de- 
viation vectors. 

The scope of this paper is to present, analyze and com- 
pare different numerical methods for the integration of 
the variational equations. In our study we consider meth- 



ods based on symplectic and non-symplectic integration 
techniques. The integration of the variational equations 
by non-symplectic methods is straightforward since one 
simply has to integrate these equations simultaneously 
with the equations of motion. This approach requires in 
general, more CPU time than schemes based on symplec- 
tic integration techniques for the same order of accuracy 
and integration time step. For this reason we focus our 
attention on methods based on symplectic schemes, ex- 
plaining in detail their theoretical foundation and apply- 
ing them to Hamiltonian systems of different numbers of 
degrees of freedom. 

The numerical solution of the variational equations ob- 
tained by the various integration schemes studied are 
used for the computation of the spectrum of the Lya- 
punov Characteristic Exponent (LCEs) and the GALIs. 
We chose to compute these two chaos indicators among 
the indices based on the evolution of deviation vec- 
tors, because the computation of the mLCE is the el- 
der and most commonly employed chaos detection tech- 
nique, while the computation of the whole spectrum of 
LCEs and GALIs requires the evolution of more than one 
deviation vector and thus is strongly influenced by inac- 
curacies of the integration procedure. We investigate the 
numerical efficiency of the different integration methods 
by comparing the CPU times they require for the com- 
putation of the LCEs and the GALIs, as well as their 
accuracy in reproducing well-known properties of these 
chaos indicators. In particular, we check whether the set 
of computed LCEs consists of pairs of values having op- 
posite signs, and if the time evolution of GALIs follows 
specific theoretically predicted laws. 

The paper is organized as follows: after introducing the 
concept of variational equations in the next section, we 
describe in sections IlIII and ITVl the LCEs and the GALIs 
respectively, which are the two chaos indicators we use 
in our study. Then, in section |V] we give the basic prop- 
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ertics of symplcctic integrators. Section [VT] is devoted to 
the detailed description of several numerical schemes for 
the integration of the variational equations of Hamilto- 
nian systems. Applications of these schemes to regular 
and chaotic orbits of systems with two or more degrees 
of freedom are presented in section IVID where also the 
efficiency of each technique is discussed. Finally, in sec- 
tion [VllTl wc summarize the results and present our con- 
clusions, while in the appendix the explicit expressions 
of the various integration methods for the Henon-Heiles 
system are given. 



II. THE VARIATIONAL EQUATIONS 

Let us consider an autonomous Hamiltonian system of 
N degrees of freedom (iVD) having a Hamiltonian func- 
tion 

H{qi,q2,. . . ,qN,Pi,P2, ■ ■ ■,Pn) = ft- = constant, (1) 

where qi and pi, i = 1,2, ... ,N are the generalized co- 
ordinates and conjugate momenta respectively. An orbit 
in the 2A^-dimcnsional phase space S of this system is 
defined by the vector 

^(i) = {qi{t),q2{t), . . . ,qN{t),pi{t),p2{t), . . . ,PN{t)), 

(2) 

with Xi = qi, Xi+N = Pi, * = 1, 2, . . . , N . The time evolu- 
tion of this orbit is governed by the Hamilton's equations 
of motion, which in matrix form are given by 



OH 
dp 



dH 
' 9q 



— J2W 



D 



(3) 



with q ^ {qi{t),q2it),. . . ,qNit)), 
{pi{t),P2it), . . . ,PN{t)), and 
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H — 



OH dH 
dqi dq2 



dH OH dH 
dqN dpi dp2 



an 

dpN 



with (^) denoting the transpose matrix. Matrix J2N has 
the following block form 



I2N 



On 

—In On 



with In being the N x N identity matrix and On being 
the N X N matrix with all its elements equal to zero. 

An initial deviation vector w{0) = 
{Sxi{0),Sx2{0), ■ . ■ ,Sx2n{0)) from an orbit x{t) evolves 
in the tangent space TsS of S according to the so-called 
variational equations 



[j2N-Tii{m)] 



A{t) -w. 



(4) 



with D|j(x(t)) being the Hessian matrix of Hamiltonian 
([1]) calculated on the reference orbit x{t), i. e. 



dxidxj 



, hi = 1,2, , 



,2N. 



Equations (|4]) are a set of linear differential equations 
with respect to w, having time dependent coefHcients 
since matrix A(t) depends on the particular reference 
orbit, which is a function of time t. 

In the present paper we consider autonomous Hamil- 
tonians of the form 



1 ^ 



(5) 



with V{q) being the potential function. The Hamilton's 
equations of motion ([3]) become 



(6) 



f 




p 






dV{q) 


.p . 




L dq \ 



while the variational equations ^ of this system take 
the form 



5q 
Sp 



= A(^) • w 



On 



In ' 




- 5q- 


On 




Sp 



6q ~ 5p 

5p = -Bl.{q{t))Sq 



(7) 



with 6q ~ 

{5pi{t),5p2{t) . . .,SpN{t)), and 
d^V{q) 



{Sqi{t),dq2it),...,dqN(t)), Sp = 
3,k^l,2,...,N. (8) 



dq-jdqk 



S(t) 



Thus, the tangent dynamics of Hamiltonian ([5|) is repre- 
sented by the time dependent Hamiltonian function 

N N 

HviSq, Sp; t)^^J2^P^ + 2T. ^HmhkSqjSqk, (9) 

1=1 " j,k 

which we call the 'tangent dynamics Hamiltonian' 
(TDH), and whose equations of motion are exactly the 
variational equations ([7]). 



III. THE LYAPUNOV CHARACTERISTIC 
EXPONENTS 

The LCEs are asymptotic measures characterizing the 
average rate of growth (or shrinking) of small perturba- 
tions to the solutions of a dynamical system. Their con- 
cept was introduced by Lyapunov when studying the sta- 
bility of non-stationary solutions of ordinary differential 
equations [l2|, and has been widely employed in study- 
ing dynamical systems since then. A detailed review of 
the theory of the LCEs, as well as of the numerical tech- 
niques developed for their computation can be found in 
i. 

The theory of LCEs was applied to characterize chaotic 
orbits by Osclcdec [13|, while the connection between 
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LCEs and exponential divergence of nearby orbits was 
given in [iJ: [3] ■ For a chaotic orbit at least one LCE is 
positive, implying exponential divergence of nearby or- 
bits, while in the case of regular orbits all LCEs are zero 
or negative. Therefore, the computation of the mLCE 
Xi is sufficient for determining the nature of an orbit, 
because xi > guarantees that the orbit is chaotic. 

The mLCE is computed as the limit for t — >■ oo of the 
quantity 

1. wmw 



■In 



t |K(0)||' 



(10) 



often called /iniie time mLCE, where ^(O), w(t) are devi- 
ation vectors from a given orbit, at times t = and t > 
respectively, and || • || denotes the norm of a vector. So, 
we have 



Xi = lim Xiit). 



(11) 



If the energy surface defined by ([T]) is compact, it has 
been shown that this limit is finite, independent of the 
choice of the metric for the phase space and converges to 
Xi for almost ah initial vectors w{0) [Ulli, [l^- Xi{t) 
tends to zero in the case of regular orbits following a 
power law [13] 

Xi(t)(xt-\ (12) 

while it tends to nonzero values in the case of chaotic 
orbits. 

An iVD Hamiltonian system has 27V (possibly non- 
distinct) LCEs, which are ordered as xi 5^ X2 > • • ■ > 
X2N- In [3 ^ theorem was formulated, which led di- 
rectly to the development of a numerical technique for 
the computation of all LCEs, based on the time evolu- 
tion of many deviation vectors, kept linearly independent 
through a Gram-Schmidt orthonormalization procedure. 
The theoretical framework, as well as the corresponding 
numerical method for the computation of all LCEs (usu- 
ally called the 'standard method'), were given in (l6l.[l7|. 
According to this method all other LCEs X2, X3 stc, 
apart from the mLCE obtained from (|lip . are computed 
as the limits for i — >■ oo of some appropriate quantities 
X2(t) , X^( t) etc., which are called the finite time LCEs 
(see 0, [l3] for more details) . We note that throughout 
the present paper, whenever we need to compute the val- 
ues of the LCEs, we apply the discrete QR-decomposition 
technique [l^. Sect. 2.10], which is a variation of the stan- 
dard method (see Sect. 6.3 of 3 for more details). 

It has been shown in [l^ that in the case of an au- 
tonomous Hamiltonian flow, the set of LCEs consists of 
pairs of values having opposite signs 



Xi 



-X2N-1+1 , i = l,2,...,N. 



(13) 



In addition, since the Hamiltonian function is an integral 
of motion, at least two LCEs vanish, i. e. 



Xn — Xn+1 — 0, 



(14) 



while the presence of any additional independent integral 
of motion leads to the vanishing of another pair of LCEs. 



IV. THE GENERALIZED ALIGNMENT INDEX 

The GALI is an efficient chaos detection technique in- 
troduced in as a generalization of a similar indica- 
tor called the smaller alignment index (SALI) [sj. The 
method has been applied successfully for the discrimina- 
tion between regular and chaotic motion, as well as for 
the detection of regular motion on low dimensional tori 
to different dynamical systems [13, [13] ■ 

The GALI of order k (Gk) is determined through the 
evolution of 2 < < 2A^ initially linearly independent 
deviation vectors WiiO), i = 1,2, . . . , k. The time evo- 
lution of each deviation vector is governed by the varia- 
tional equations ([7]) . Each evolved deviation vector Wi (t) 
is normalized from time to time, having its norm equal to 
1, in order to avoid overflow problems, but its direction 
is left intact. Then, according to [3], Gk is defined to be 
the volume of the fc-parallelogram having as edges the k 
unitary deviation vectors Wi{t), i = 1, 2, . . . , fc. This vol- 
ume is equal to the norm of the wedge product of these 
vectors, and Gfc is given by 

Gk{t) = WMi) A Mi) A • • • A Wk{t)\\- (15) 

From this definition it is evident that if at least two of the 
deviation vectors become linearly dependent, the wedge 
product in (|15|) becomes zero and the Gk vanishes. 

Expanding the wedge product (jlSp into a sum of deter- 
minants and studying the asymptotic behavior of those 
who vary the slowest in time, it is possible to show ana- 
lytically the following [3]: in the case of a chaotic orbit 
all deviation vectors tend to become linearly dependent, 
aligning in the direction defined by the mLCE and Gk 
tends to zero exponentially following the law 



Gk{t) (X e 



-[(cri-(T2) + (cri-(T3)H V{<Ti-ak)]t 



(16) 



where cti , . . . , tTfe are approximations of the first k largest 
Lyapunov exponents. On the other hand, in the case of 
regular motion on an A'^-dimensional torus, all deviation 
vectors tend to fall on the iV-dimensional tangent space 
of this torus. Thus, if we start with k < N general devia- 
tion vectors they will remain linearly independent on the 
A^-dimensional tangent space of the torus, since there is 
no particular reason for them to become aligned. As a 
consequence Gk is different from zero and remains prac- 
tically constant for k < N. On the other hand, Gk tends 
to zero for k > N, since some deviation vectors will even- 
tually become linearly dependent, following a particular 
power law which depends on the dimensionality N of the 
torus and the number fc of deviation vectors. The behav- 
ior of Gk for regular orbits lying on iV-dimensional tori 
is given by 

„ , , f constant if 2 < fc < 
^'^('^"l^^ ffA^<fc<2iV (17) 

If the regular orbit lies on a low dimensional torus, i. e. 
an s-dimensional torus with 2 < s < N then Gk remains 
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practically constant and different from zero for k < s and 
tends to zero for k > s following particular power laws 
(see for more details). 

In order to compute the value of Gk we consider the 
2N X k matrix W(t) having as columns the coordi- 
nates Wji{t) of the unitary deviation vectors Wi(i), i = 
1, 2, . . . , /c, j = 1, 2, . . . , 2N, with respect to the usual or- 
thonormal basis ei = (1, 0, 0, ... , 0), 62 = (0,1,0,..., 0), 
e2N = (0, 0, 0, ... , 1) of the 2iV-dimensional tangent 
space TxS and perform the Singular Value Decomposi- 
tion (SVD) of this matrix. Then, as it was shown in 
(lOj . Gk is equal to the product of the singular values z^, 
« = 1, 2, . . . , of matrix W(t), i. e. 



Gkit)^llz.,it) 



SYMPLECTIC INTEGRATORS 



Let us discuss in some detail how we can integrate the 
equations of motion of a general Hamiltonian ([T]) by 
a symplectic integration scheme, focusing our attention 
on a particular family of integrators presented in (2l| . 
Defining the Poisson bracket of functions f{q,p), gi<l,P) 
by M - 



of products of e'^''^^'^ and e'^''^^^. In this context sev- 
eral symplectic integrators of different orders have been 
developed by various researchers [H, [23] . 

In [2l| the families of SBAB (and SABA) symplectic 
integrators, which involve only forward (positive) inte- 
gration steps were introduced. These integrators were 
adapted for the integration of perturbed Hamiltonians of 
the form H ~ A + eB, where both A and B are inte- 
grable and e is a small parameter. A particular integra- 
tor SBAB„ {SBn), or SABA„ {Sau), involves n steps, 
i. e. n applications of products of e'^''^^^ and e'**'^^'^, 
and is of order C'(r^"e -I- r^e^) with respect to the in- 
tegration step T. This means that by using these in- 
tegrators, we are actually approximating the dynamical 
behavior of the real Hamiltonian ^ -I- ei? by a Hamilto- 
(18) nian H* = A + eB + 0{T^'^e + r^e^), i. e. we introduce 
an error term of the order T^"e + r'^e^. 

The accuracy of the Ssn (SAn) integrator can be im- 
proved when the commutator term C = {B , {B , A}} [2^ 
leads to an integrable system, as in the common situation 
of A being quadratic in momenta p and B depending only 
on positions q. In this case, two corrector terms of small 
backward (negative) steps can be added to the integrator 
Sbu 



N 



^ V 9qi dpi dpi dqi 



(19) 



the Hamilton's equations of motion ^ take the form 



^ = {x, H} = Lhx, 



(20) 



where Lh is the differential operator defined by L-^f ~ 
{fiX}- The solution of Eq. for initial conditions 

x(0) = xq, is formally written as 



j-ri 



(21) 



n>0 



Let us assume that the Hamiltonian function H can 
be split into two integrable parts as H = A + B. A sym- 
plectic scheme for integrating equations (PU)) from time 
t to time t + T consists of approximating, in a symplec- 
tic way, the operator e^^" = q'^C^a+Lb) integra- 
tor of j steps involving products of operators ef-^"^^-^ and 
gdiTLs^ z = 1, 2, . . . , J, which are exact integrations over 
times CiT and diT of the integrable Hamiltonians A and 
B. The constants Ci, d^, which in general can be posi- 
tive or negative, are chosen to increase the order of the 
remainder of this approximation. So e"^^-^, e^^^ are ac- 
tually symplectic maps acting on the coordinate vector x. 
Therefore the integration of equations pO| over one time 
step r, which evolves the initial coordinate vector x{t) to 
its final state x{t -f r), is represented by the action on 
x{t) of a symplectic map S produced by the composition 



SI 



Lc 



(22) 



A similar expression is valid also for Sau- The value 
of constant g is chosen in order to eliminate the r^e^ 
dependence of the remainder which becomes of order 
0(r2"e-Hr4e2). The SBAB (SABA) integrators have al- 
ready proved to be very efficient for the numerical study 
of different dynamical systems [21I, [2^, . We note that 
several authors have used commutators for im pro ving the 
efficiency of symplectic integrators (e. g. [28l. |29|). 

Setting e = 1 we can apply the SBAB (SABA) integra- 
tion schemes for the integration of Hamiltonian ([5]), since 
this Hamiltonian can be written a& H = A + B, with 



1 ^ 

A{p) = -Y,pI B{q)^V{q), 



(23) 



being both integrable. The maps e^^^ , e^^^ , which prop- 
agate the set of initial conditions (9,p) at time to their 
final values {q' ^p') at time t -\- t, for the Hamiltonian 
functions A{p) and B{q) ((23|) are 



and 



rL^ .19 = q + pT 

p = p 



q = q 

■ ^ ^, ^ dV(q) 
P = P K^''' 



(24) 



(25) 



respectively. For Hamiltonian ([5]) the corrector term is 
given by 

C = (B.iB,A}}^t(?m)\ (26, 
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which is a function of only the coordinates q and thus 
easily integrated as 



Integration of the tangent dynamics 
Hamiltonian 



. dC{q) 
P 

oq 



(27) 



In Appendix lA II we give the explicit formulas of equa- 
tions ([24]), ([25]) and ([271) for the Henon-Heiles system 

(El. 



VI. NUMERICAL INTEGRATION OF 
VARIATIONAL EQUATIONS 



In this section we present several numerical schemes 
for the integration of the variational equations, consider- 
ing both non-symplectic techniques and methods based 
on symplectic integrators. The latter schemes are quite 
general and any symplectic integrator can be used for 
their implementation. In our study we consider an effi- 
cient fourth order symplectic integrator, the S%2 [2ll.[28l|. 
which has an extra degree of complexity with respect 
to integrators composed of products of maps e"^^-*, and 
e^^^ , since it requires the application of the corrector 
term C 



A. Non-symplectic schemes 

In order to follow the evolution of a deviation vector, 
the variational equations ([7]) have to be integrated simul- 
taneously with the Hamilton's equations of motion ([SJ, 
since matrix TDyit) depends on the particular reference 
orbit x{t), which is a solution of equations ([6]). Any non- 
symplcctic numerical integration algorithm can be used 
for the integration of the whole set of equations ^ and 
©• 

In our study wc use the DOP853 integration method 
which has been proven to be very efficient. The DOP853 
integrator [soj is an explicit non-symplectic Runge-Kutta 
integration scheme of order 8, based on the method of 
Dormand and Price (see [l|, Sect. II. 5]). Two free param- 
eters, T and i5, are used to control the numerical perfor- 
mance of the method. The first one defines the time span 
between two successive outputs of the computed solution. 
After each step of length r the values of LCEs (GALIs) 
are computed and the deviation vectors are orthonormal- 
ized (normalized). For the duration of each step r, the 
integrator adjusts its own internal time step, so that the 
local one-step error is kept smaller than the user-defined 
threshold value 6. For DOP853 the estimation of this lo- 
cal error and the step size control is based on embedded 
formulas of orders 5 and 3. 



Another approach to compute the evolution of devia- 
tion vectors is to initially integrate the Hamilton's equa- 
tions of motion ([6|) , in order to obtain the time evolution 
of the reference orbit x(t), and then to use this numeri- 
cally known solution for solving the equations of motion 
of the TDH ([5]), which are actually the variational equa- 
tions dll). 

In practice one numerically solves the Hamilton's equa- 
tions of motion ^ by any (symplectic or non-symplectic) 
integration scheme to obtain the values x{ti) at ti = i At, 
i = 0,1,2,..., where At is the integration time step of 
these orbits. Of course, the accuracy of the particular 
numerical scheme used for the construction of the time 
series x{ti) will affect the quality of the numerical solu- 
tion of the variational equations, regardless of the numer- 
ical scheme used for solving them. Having computed the 
values x{ti) different methods can be applied for approx- 
imating the solution of the variational equations, which 
will be discussed in the following sections. 



1. TDH with piecewise constant coefficients 

One method is to approximate the actual time depen- 
dent TDH ([9]) by a Hamiltonian with piecewise constant 
coefficients. This means to assume that the coefficients 
D^(g(t))jfe j, fc = 1,2, ...,A^ of iJv ^ are constants 
equal to 'Dy{q{ti))jk for the time interval [ti,ti + /S.t). 
These constants are determined by the values of the or- 
bit's coordinates and are known, since we know the time 
series x{ti) = {q{ti),p{ti)). Thus, for each time interval 
[ti, ti + At) wc end up with a quadratic form Hamiltonian 
function Hv{Sq, 6p;ti), whose equations of motion form 
a linear system of differential equations with constant 
coefficients. 

The Hamiltonian Hv{5q,5p;ti) can be integrated by 
any symplectic or non-symplectic integration scheme, or 
can be explicitly solved by performing a canonical trans- 
formation to new variables Q, P, so that the trans- 
formed Hamiltonian HyQp becomes a sum of uncoupled 
ID Hamiltonians, whose equations of motion can be in- 
tegrated immediately. To this end, let \k be the eigen- 
values and k = 1,2, . . . , N the unitary eigenvectors 
of the constant matrix 'Dy{q{ti)). Then matrix T, hav- 
ing as columns the eigenvectors v^, defines a canonical 
change of variables q = TQ, p = TP, which gives Hy 
the diagonal form 



H 



VQP 



^ 1 



(28) 



The equations of motion of Hyqp are then easily solved. 

In our study we use the same symplectic integrator 
(S'ljg) both for obtaining the time series x{ti) and for in- 
tegrating the quadratic form Hamiltonian Hy(6q, Sp;ti) 
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in the time interval [ti,ti + At). We name this approach 
the TDHcc method (cc: constant coefficients). An alter- 
native approach is to compute the exact solution of the 
equations of motion of Hv{Sq, Sp;ti) (whose piecewise 
constant coefficients are obtained by the symplectic inte- 
gration of the orbit using the scheme) by transform- 
ing it to a system of N uncoupled harmonic oscillators 
through the canonical transformation induced by matrix 
T. This approach is called the TDHes method (es: exact 
solution) . 

In general, the transformation matrix T is determined 
for each time interval [ti, ti + At) by solving numerically 
the eigenvalue problem 



The dynamics of the ND TDH Hy & is equivalent to 
that of the (TV + 1)D Hamiltonian 



~ ^ ^ 1 ^ . 



1 ^ 



(32) 



This Hamiltonian can be easily integrated by any sym- 
plectic integration scheme, since it can be split into two 
integrable parts 



T>liq{t,))v = Xv, 



(29) 



a procedure which could become computationally very 
time consuming, especially for systems with many de- 
grees of freedom. On the other hand, in some simple low 
dimensional cases, like for example the Henon-Heiles sys- 
tem (|54p . the transformation matrix T can be determined 
analytically (see Appendix IA2ap . 



2. Integration of the TDH in an extended phase space 

Instead of approximating Hy (H]) by a quadratic 
form having constant coefhcients for each time interval 
[ti,ti + At), we can explicitly treat Hy as a time de- 
pendent Hamiltonian. This time dependency is due to 
the fact that the coefficients of Hy are functions of the 
orbit's coordinates q{t). Like in the previous approach, 
we consider the time series q{ti) to be known from the 
numerical integration of the Hamilton's equations ([6|). 

The iVD time dependent Hamiltonian Hy can be 
transformed to a time independent Hamiltonian Hy with 
an extra degree of freedom by considerin g th e time t as an 
additional coordinate (see for example [3l|, Sect. 1.2b]). 
For this purpose, we add to the Hamilton's equations of 
motion of Hy the equations 



t = 1 , Hv 



dt 



(30) 



Then we set t and ~Hy as an additional coordinate and 
momentum respectively, i. e. ^^at+i = t, 5pjv+i = ^Hy, 
and define the new Hamiltonian 



Hy iCv) ^ Hy{Sq,Sp:t) + SpN^ 



(31) 



where ^ = {5q,t) and ff ^ (6p, —Hy) are respectively the 
new coordinates and momenta. The flow in the {2N + 
2)-dimensional extended phase space of the {N + 1)D 
Hamiltonian Hy is parameterized by a 'new' time such 
that i(^) — C, which does not appear explicitly in the 
functional form of Hy ([51]) . The set of equations ([7]) and 
are the Hamilton's equations of motion of Hy. 



1 ^ 

A{Sp,SpN+l) = - ^ (5p,? -f (5pAr+i 



N 



(33) 



B{bq,t) = \ Y,^\m))j}^HM. 



The maps e'^^-i, e'^^s. which propagate the set of initial 
conditions {5q, t, 6p, SpN^i) at time t, to their final values 
{Sq ,t' ,5p , Sp'j^j^i) at time t + r are 



„rLi 



Sq = Sq + SpT 

tL, . ^t'^^ = t + T 

Sp = 5p 

= SpN+l 



6q = 6q 
t' = t 

dB{5q,t) 

op = op T 



(34) 



^Pn+i = ^P 



N+l 



dSq_^ 
dB{Sq,t) 

dt ' 



(35) 



The corrector term of the SBAB and SABA integration 
schemes 



N 



,^1 V dSq^ ^ 



(36) 



is a function of only the coordinates ^ = (6q, t) and thus 
easily integrated 



e^^e : < 



Sq = Sq 

ip Jsp-'M^r 



dSq^ 

dC{5q,t) 

dt ' 



(37) 



The explicit expressions of these maps for the Henon- 
Heiles system ((M)) arc given in Appendix lA 2 bl 
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From equations (|34)) . ((35|) and ((37)) we see that time 
t is changed only by the act of operator e"^^-?. On the 
other hand, operators e'^^s and e'^^c require the knowl- 
edge of positions q at specific times for the evaluation of 
the partial derivatives of B and C. We also note that 
for all these operators the last equation for SpN+i can 
be neglected, since the knowledge of its value does not 
influence the evolution of the other quantities, and con- 
sequently the solution of the variational equations ([7]). 

Since the coordinates of the orbit q are known only at 
specific times ti = iAt, i = 0,1, . . ., one is restricted to 
use integration schemes that require the knowledge of q 
at exactly these times. Such a scheme is, for example, 
the Sbi integrator 



Sbi = e,(-/2)^Be^i-4e(-/2)i^ 



(38) 



(which is practically the well-known Stormer/Verlet or 
leap-frog method) with r = At. The right opera- 
tor e^'^/^^^B which acts first, requires the knowledge of 
q{ti), while the left operator e^'^/^'^s needs the values of 
q{ti + t) = q{ti+i), because the time value has changed 
from ti to ti + T by e"^^-? . Note that the 5*^1 integrator 



Sai 



,(r/2)L^ rig (t/2)L^ 



(39) 



requires the knowledge of q{ti + t/2) for the application 
of e'^^B. This second order integration scheme could be 
used with t = 2At, leading in general to a less accurate 
algorithm compared to Sbi (PS]) . which is also a second 
order integrator but uses a smaller time step r = At. 
For r = 2At is in general more efficient to apply the 
integration scheme 



Sb2 = e(^/^)^Be(^/2)L^g(2r/3)Lgg(r/2)L^g(r/6)L^ 



(40) 



which was initially derived in |24| . This integrator needs 
the known values q{ti), q{ti +t/2) = q{ti + At) = q{ti+i) 
and q{t, + r) = q{t, + 2 At) = q{t,+2)- 

The above integration schemes can also be combined 
with a corrector step, since e^^e ([57)) does not change 
the time values, and acts before and after the main body 
of the integrator (see equation ([H])), when t has values 
for which we know the coordinates q. We refer to this 
technique as the TDHeps method (eps: extended phase 
space). For the numerical applications of the TDHeps 
method (presented in Sect. I VII)) we use the fourth order 
integrator Sb2 both for the integration of the variational 
equations and for the computation of the orbit. 

Higher order SBAB or SABA integrators cannot be 
used in this framework, because they require the knowl- 
edge of q at non equidistant time values, different from 
ti. In order to apply such schemes one could initially 
compute the solution of equations (|6]) also at these spe- 
cific times (e. g. by interpolation), but this would lead 
to a cumbersome, complex, time consuming, and conse- 
quently inefficient scheme. 



C. The tangent map (TM) method 

The set of equations (|6]) and ([7)) can be considered as 
a unified set of differential equations 



P 



P 



dV{q) 



Sq = Sp 

Sp^-Dl{q)5q } 



du 
H 



JHVU, 



(41) 



where u ~ {q,p, Sq, 6p) is a vector formed by the phase 
space vector x = ((f, p) and the deviation vector w ~ 
{6q, Sp), and Lhv is the differential operator of the whole 
system. In analogy to equation ((^T)) . the solution of sys- 
tem ()4T)) for an initial condition u{0) can be formally 
written as u{t) = e*^"^M(0). We describe now how sym- 
plcctic integrators can be used to obtain this solution. 

First of all, let us note that equations ()^T)) cannot be 
considered as the Hamilton's equations of motion of some 
generalized Hamiltonian function. If such a Hamiltonian 
existed, and could be split into two integrable parts, any 
symplectic integrator could be used for finding the solu- 
tion of system ()4T)) . Since this is not the case, we follow a 
different approach to achieve this goal. In section |V] the 
integration of the equations of motion of Hamiltonian ()5)) 
over one integration time step r was split into steps over 
appropriate time intervals c^r, diT, where the dynamics 
was determined either by Hamiltonian A{p) or B(q) ()23p . 
During these intermediate steps the tangent dynamics of 
the system is governed by the variational equations 



for A{p), and by 



Sq = Sp 
Sp^O 



Sq = Q 

s'p=-Dl.(q)Sli 



(42) 



(43) 



for B{q). Therefore, for each intermediate step of the 
symplectic integration scheme the dynamics of the phase 
and the tangent space is governed by the set of equations 




Lavu, 



(44) 



and 



q = 



P 



dV{q) 



dq 

5g = 

5'i>=-T>l{q)5q } 



du ^ 
-dt = 



(45) 
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for Hamiltonians A{p) and B{q) (|23|) respectively, with 
Lav Slid Lgy being the corresponding differential oper- 
ators. 

These sets of equations are immediately solved, leading 
to maps 




(46) 



q' 


= q 


P' 


= V- 


5q 




/ 

5p 


= 5p 



dV{q) 



(47) 



Obviously the first two equations of maps e'^^^^' , e"^^^^' 
are exactly maps e^^^ ([M]) and e'^^^ (P5|) . respectively. 

Thus, any symplectic integration scheme used 
to solve the Hamilton's equations of motion ([6|), 
v^rhich involves the successive application of maps 
e^L* ([24)) . e^^Ei ([25])^ can also be used for the si- 
multaneous integration of the variational equa- 
tions ([7]), i. e. for solving the set of equations 
(|4ip. by replacing maps e^^^^, e^'^^' with maps 
e^LAv ([46]) and e^^^v respectively. This state- 

ment is a specific application of a more general result 
which is stated for example in [2l|: Symplectic inte- 
gration schemes can be applied to first order differen- 
tial systems X = LX that can be written in the form 
X = (La + Lb)X, where i. La and Lb are differen- 
tial operators for which the two systems X = LaX and 
X = LbX are integrable. The system of differential 
equations u = Lhvu (j4ip belongs to this category since 
it can be split into the integrable systems ii = Lavu p4)) 
and u ~ Lgyu (j45p . 

Let us discuss this splitting in more detail. The system 
(|^T]) can be written as 



Q = V 



(48) 



with Q ^^{q,Sq) = {qi,q2, . . . ,qN,5qi,5q2, . ■ . ,SqN), 
V = {p,Sp) = {pi,p2, . ■ ■ ,PN,Spi,Sp2, . . . ,6pn), and 
J^{Q) being a vector with coordinates 



given by 




d 



AV 



^BV 



)u{Q,r). 



(50) 

The solution of Eq. ([50]) for a time step r can be formally 
written as 

U{t + t)= e^^^^'-+^^''^U{t). (51) 

The decomposition of e'^(^'*^'+^^^) into products of oper- 
ators e'^'"^^ , e^^^^ by any symplectic integration scheme 
gives rise to an exponential-splitting algorithm for the 
integration of system (HI]), which would be symplectic if 
Eqs. (|4T|) were the equations of motion of a Hamiltonian 
function (which are not, as we have already discussed). 

In our study we consider symplectic integrators that re- 
quire the application of corrector terms. When the 5*1;^^ 
i^An) integrators are used, map e^^'^ (|27)) acts for some 
intermediate steps of the algorithm. Formally one can 
consider that for these steps the phase space dynamics 
is governed by the Hamilton's equations of motion of the 
Hamiltonian function C{q) ([26]) (whose solution is given 
by map e^^'^ p7p ). Consequently, the tangent space dy- 
namics is described for these time steps by the variational 
equations of Hamiltonian C{q). So the evolution of the 
general vector u is given by 



dC{q) 



dq 



P = 

= 

5p=-'Dl{q)5q J 



du 
It 



^cvu, 



(52) 



where D^(g)jfe = d'^C[q) / dqjdqk- We easily see that the 
solution of these equations is given by the map 



„tLcv 





= q 


P' 


= p- 


5q 


= Sq 


. 5p 


= 6p - 



dCjq) 
dq 



nUq)SqT 



(53) 



dqi 



for 1< i < iV, 



which, of course, is an extension of map e^^^ ([27]). So 
the use of the corrector term with the Sen (Sah) 
integrator for the integration of system (|4ip re- 
quires the additional substitution of map e^^"^ 
([27]) by the extended map e^^«=v ([53]) . 

We call the above-described procedure for the simulta- 
neous integration of the Hamilton's equations of motion 
([6]) and the variational equations ([7]), the tangent map 
(TM) method. The explicit expressions of the extended 
maps e^-^-*^ ([H]), e^-^^^ ^ and e^-^^^ ([S3]) for the 
Then the dynamics of any general variable U{Q,V) is Hcnon-Heiles system ([M]) are given in Appendix lA 2 el 



j^^X^Sq^ for 7V<^<27V 



(49) 



fc=i 



dqidqk 
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VII. NUMERICAL APPLICATIONS 

In order to study the efficiency of the different schemes 
for the integration of the variational equations, we apply 
them to some simple Hamiltonian systems of different 
numbers of degrees of freedom. In particular we consider 

a) the well-known 2D Henon-Heiles system [s^l described 
by the Hamiltonian 

H2 = \{pI +pI) + lix' + y') + x'y ly^ (54) 

b) the 3D Hamiltonian system 



0.4 - 



■i 2^ 



2\ , 2 
-Px) + —(y fyj 

(55) 

studied in 0, [13, [s^l, and c) the famous Fermi-Pasta- 
Ulam (FPU) /3- lattice model [SJI , which describes a chain 
of N particles with nearest neighbor interaction, for the 
particular case of A'' = 8 studied in The 8D Hamil- 
tonian of this system is 



^ 2 

i=l 



8 

E 

i=0 



(56) 

We consider some typical regular and chaotic orbits of 
these systems and investigate the efficiency of the vari- 
ous numerical techniques by checking how well their out- 
comes verify the following theoretically known properties 
of the LCEs and the GALIs: 

• The finite time mLCE Xi (t) should eventually tend 
to zero in the case of regular orbits following the 
power law given in (jI2p . 

• According to Eq. p^ . the LCEs are grouped in 
pairs of values having opposite signs, and conse- 
quently their sum vanishes. Therefore the same re- 
lation should be also satisfied by the limiting values 
of the corresponding finite time LCEs i. e. 

\im {X,{t)+X2N-^+l{t))=0 , i = l,2,...,Af. (57) 



• According to Eq. (|T4)) at least two LCEs vanish and 
therefore A'jv(i) and A'jv+i(i) should tend to zero. 

• The GALIs follow the laws ([TH) and for chaotic 
and regular orbits respectively. 



A. The 2D Henon-Heiles system 

We implement first the various numerical schemes pre- 
sented in Sect. IVII for the integration of the variational 
equations of regular and chaotic orbits of the 2D Henon- 
Heiles system ([M)) . The explicit expressions of all these 
schemes are presented in detail in Appendix lA 21 The 
orbits of the Henon-Heiles system have four LCEs xi > 



-0.4 - 



- * 




FIG. 1: The PSS defined by s = 0, > 0, for the Henon- 
Heiles system ([54]) with H2 = 0.125. The regular orbit Rl 
corresponds to the five closed black curves around the right 
large island of stability, while the chaotic orbit CI is repre- 
sented by the black dots scattered over the PSS. In order to 
get a clear picture of the structure of the whole PSS, other 
orbits of the system are plotted in gray. 



X2 > X3 > X4, with X2 = X3 = and xi = -Xi > 0. 
A simple qualitative way of studying the dynamics of a 
Hamiltonian system is to plot the successive intersections 
of its orbits with a Poincare surface of section (PSS) (see 
for example Sect. 1.2b of [HI). In 2D systems like ((54| . 
the PSS is a two dimensional plane which allows the clear 
visualization of the dynamics. 

In our study we keep the value of the Hamiltonian 
fixed at H2 = 0.125. Initially, we consider two repre- 
sentative orbits of the system: the regular orbit Rl with 
initial conditions x = 0, px ^ 0.2334, y = 0.558, Py = 0, 
and the chaotic orbit CI with initial conditions x = 0, 
Px « 0.4208, y = -0.25, Py = 0. In Fig. [1] we plot the 
intersection points of these two orbits with the PSS de- 
fined by a; = 0, > 0. The points of the regular orbit 
lie on a torus and form five smooth closed curves (the 
so-called stability islands) on the PSS, while the points 
of the chaotic orbit appear randomly scattered. 

First, we use the DOP853 non-symplectic scheme to 
integrate the set of differential equations composed from 
the Hamilton's equations of motion (|A1[) and the varia- 
tional equations (|A2|) . In our computations we set the 
integration time step r = 0.05 and the threshold param- 
eter S = 10~^, unless otherwise stated. 

We also implement the TDHcc, the TDHes and the 
TDHeps methods. For these methods we initially inte- 
grate equations (|Aip by the scheme. In this way 
we obtain the coordinates of the orbit at times tt = iAt, 
i = 0,1,2,..., with At being the constant integration 
step. Then we assume the TDH (|A3[) to have constant 
coefficients in each time interval [ti, ti + At) and either we 
integrate in this interval its equations of motion by the 
integrator (TDHcc method), or we compute the ex- 
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act solution of these equations by performing the canon- 
ical transformation induced by matrix T of Eq. (|A10|) 
(TDHes method). Alternatively, we use the scheme 
for integrating the equations of motion of the 3D Hamil- 
tonian Hvh (|A13p in the time interval [ti,ti + 2 At), by 
applying Eqs. (jAlSj) . (|A16P and (|A18|) with time step 
T = 2 At (TDHcps method). Finally we implement the 
TM method using the integrator, which requires the 
application of maps (|A21|) . (|A22p and (|A23|) . 

As a final remark we note that in all the above- 
described schemes after each time step r the LCEs 
(GALIs) are computed and the deviation vectors are or- 
thonormalized (normalized) having norm equal to 1. 



1. Regular orbits 

Results concerning the LCEs of the regular orbit Rl 
are shown in Fig. [5] In particular, the time evolution of 
the finite time LCEs Xi and X2 is given in the upper 
panels, while in the lower panels the evolution of quanti- 
ties \Xi +X4\, -h X3I is plotted. 

In Table U information on the computation of the whole 
spectrum of LCEs of the Rl orbit up to < = 10^ is re- 
ported. The relative energy error, which could be con- 
sidered as an indicator of the goodness of the integra- 
tion procedure of orbit Rl, increases with time for the 
DOP853 method, while it fluctuates around a constant 
value for all other methods. The values of this error and 
of Xi at the end of the integration arc reported in the ta- 
ble. The CPU time needed on an ordinary personal com- 
puter by each method for the integration of the equations 
of motion and the variational equations, as well as for the 
computation of the spectrum of LCEs is also given. 

The results of Fig.|l]show that the DOP853 (Fig.UJa)) 
and the TM method (Fig. [2Kc)) have the best perfor- 
mance in evaluating the mLCE, because Xi tends to zero 
until the end time t = 10® of the integration, following 
a t~^ law. The good behavior of the DOP853 and the 
TM methods is due to the fact that the first technique 
is used for the integration of the actual set of Eqs. (|A1[) 
and (|A2p which govern the dynamics of the orbit and the 
deviation vector, while the second method approximates 
very accurately the dynamics of the system by the re- 
peated application of a symplectic map, and the tangent 
dynamics by the act of the corresponding tangent map. 

For the TDHcc (Fig.[2tb)), the TDHes (Fig.I^c)) and 
the TDHcps (Fig. [D^d)) methods Xi initially decreases 
too as Xi cx but later its value deviates from the ap- 
proximate t^^ law and tends to a constant (different for 
each method) nonzero value. Among these techniques the 
TDHeps method has the best performance, because the 
computed Xi levels off to smaller values than in the cases 
of TDHcc and TDHes methods, being Xi « 2.3 x 10"^ 
at t = 10®. Nevertheless, from the results of Figs. [IJb)- 
(d) one would wrongly characterize the regular orbit Rl 
as chaotic. Concerning the TDHcc and TDHes meth- 
ods, the main reason for this discrepancy is that these 



methods approximate the tangent dynamics by consid- 
ering constants the actual time dependent coefficients of 
Hamiltonian Hyn (|A3p . for the duration of each integra- 
tion time step. The equations of motion of Hyn with 
constant coefficients are solved explicitly by the TDHes 
method; while their solution is approximated by the ap- 
plication of the TDHcc scheme. For the used time step 
T = 0.05, both methods give practically the same Xi at 
t ~ 10®. For smaller time steps the final values of Xi 
obtained by both techniques are closer to the theoretical 
value Xi = 0. On the other hand, since the TDHeps 
method takes into account the time dependent nature of 
the coefficients of Hvh, it succeeds in obtaining a better 
estimation of the mLCE compared to the TDHcc and the 
TDHes methods. 

The computed values of the second largest LCE (x2 = 
0) have similar characteristics with the results for the 
mLCE. Again, the finite time LCE X2 computed by 
the DOP853 integrator (Fig. ^a.)) and the TM method 
(Fig. [U^e)) tends to zero until the end of the integration 
time. On the other hand, the X2 computed by the TD- 
Hcc (Fig.lHb)), the TDHes (Fig.[li;c)) and the TDHeps 
(Fig-Eld)) methods does not tend to zero, but levels off 
to positive values which are always smaller than the level 
off values of Xi . Again the TDHeps approach is more ac- 
curate, because the final value X2 ~ 9.3 x 10~^ att = 10® 
obtained by this method is slightly smaller than the ones 
found by the TDHcc and the TDHes methods, and thus 
closer to the real X2 = value. 

The ability of the DOP853 and the TM methods to 
evaluate quite accurately the LCEs of the regular orbit 
Rl is also shown by the tendency of quantities -I-X4I, 
\X2 + X^l to become zero (Fig. EJf) and (j)). Actually 
these quantities attain, for both methods, very small val- 
ues < 10^'' aX t = 10®. But when these quantities are 
computed by the other three techniques they do not be- 
come zero as they theoretically should do, but level off to 
small positive values (Figs. [2I^g)-(i)). Again the TDHeps 
method exhibits a better performance since the level off 
values are smaller than the ones obtained by the TDHcc 
and the TDHes methods. 

Looking in Table |T] at the CPU times needed for the 
computation of the whole spectrum of LCEs, one sees 
that the non-symplectic method is the most expensive 
one. Amongst the remaining approaches the TM method 
is the fastest, due to the fact that the whole set of equa- 
tions for the evolution of both the orbit and the deviation 
vector are integrated together. The TDHcc and TDHes 
methods require more CPU time than the TM method, 
because for each integration time step the evolutions of 
the orbit and the deviation vectors are not performed si- 
multaneously. First the orbit is evolved. Its coordinates 
define the coefficients of Hvh (|A3p . which are consid- 
ered to be constant for the duration of the time step. 
Then, the deviation vectors are advanced for this particu- 
lar Hamiltonian function for one time step. The TDHeps 
method needs even more CPU time mainly because the 
orbit is integrated with half time step (At = t/2) with 
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DOP853 TDHcc TDHes TDHeps JM 




1234567 1234567 1234567 1234567 1234567 
loQiot loQiot log^ot logiot loQiot 

FIG. 2: The time evolution of Xi{t) (black curves), X2{t) (gray curves) [upper panels] and \Xi{t) + X4{t)\ (black curves), 
(t) + ^3(^)1 (gray curves) [lower panels] in log- log plots for the regular orbit Rl of the Henon-Heiles system ((54]). The 
variational equations are integrated by the DOP853 integrator ((a) and (f)), and by the TDHcc ((b) and (g)), the TDHes 
((c) and (h)), the TDHeps ((d) and (i)) and the TM ((e) and (j)) method. Dashed lines in panels (a) and (e) correspond to 
functions proportional to t~^. The step size is r = 0.05 for all methods. For the DOP853 method the parameter 5 = 10"' is 
used. 
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TABLE I: Information for the computation of the whole spectrum of LCEs for the regular orbit Rl of the Henon-Heiles system 
(|54|l . up to t = 10*. The non-symplectic DOP853 algorithm and the symplectic Sb2 integrator are used. In the latter case the 
Sb2 scheme is used for the evolution of the orbit, while different approaches are applied for the integration of the variational 
equations. Step size r is the time between two successive evaluations of the LCEs. For the TDHcc, the TDHes and the TM 
methods, t coincides with the integration time step At of the orbit, while for the TDHeps method r — 2At. In the case of the 
DOP853 algorithm the integration over time r is performed with a variable integration step, so that the local one-step error is 
kept smaller than 5. The relative energy error and the estimated value Xi of the mLCE at t — 10** are given. The required 
CPU time for the implementation of each method on an ordinary personal computer (AMD Athlon IGHz) is given in the last 
column. The first 5 cases (above the horizontal line) are the ones presented in Fig. (2] 



respect to the other methods. we see that the energy error for the DOP853 method at 

The first five rows of Table U contain information for t = 10*, is smaller than the error of the S%2 integrator 
the particular cases shown in Fig. O From these data used by the other methods. As it is also shown in Fig. [2] 
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FIG. 3: The time evolution of Xi{t) in log- log plots for the 
regular orbit Rl of the Henon-Heiles system (|54p for (a) the 
DOP853 (with 5 = 10"^), and (b) the TDHeps methods, 
when different step sizes r are used. In (a) the curves for 
r = 0.05 and r = 0.1 practically overlap. 



the values of Xi obtained by the DOP853 and the TM 
methods arc close to each other, despite the fact that 
the DOP853 method integrates orbit Rl with a better 
accuracy. Of major practical importance is the fact that 
the DOP853 method needs almost two times more CPU 
time than the TM method in order to compute the four 
LCEs up to i = 10®. Increasing the integration step size 
of DOP853 to r = 0.1 (Fig.E^a)) stiU permits the com- 
putation of the same Xi value at t = 10®, but with a 
larger error in the conservation of H2 ■ The Xi computed 
by the DOP853 method for even larger step sizes, like 
T = 0.2 and r = 0.5, starts after some time to exhibit 
deviations from the Xi cx law (Fig. EJa)), leading 
to somewhat larger final values {Xi « 2.4 x 10^*" for 
T = 0.2 and Xi « 1.1 x lO"** for t = 0.5) with respect 
to the Xi « 1.6 X 10~^ value found for smaller r. From 
our numerical experiments we see that the required CPU 
time for the DOP853 method, as well as the relative er- 
ror of the computed energy H2 mainly depend on the 
integration time step r and not on the threshold param- 
eter 6. In particular, for r < 0.2 the value of S does not 
practically influence the required CPU time. For larger 
values of r (for which nevertheless the obtained results 
are not very accurate) the CPU time is increased and the 
accuracy is improved when S is decreased. On the other 
hand, the TM method succeeds even for r = 0.5 to com- 
pute very fast the correct small final value of Xi < 10" 
This method keeps also the relative energy error at an 
acceptably low level, which is not the case any more for 
the DOP853 method with the same time step. Besides 
the computation speed, this is an additional advantage 
of the TM method over the DOP853 scheme. 

It is worth noting that, although the DOP853 algo- 
rithm is an integration scheme of higher order than the 
S'|j2 symplectic integrator used in the TM method, it 
shows worse characteristics than the TM method, not 
only for large t, but also when we compare implementa- 
tions of the two algorithms that require almost the same 
CPU time. For example, the DOP853 method for t = 0.2 
and S = 10"^" (or even S ~ 10~^) has a final relative en- 
ergy error which is larger by 2 orders of magnitude with 



respect to the error of the TM method for r = 0.1 (which 
requires almost the same CPU time sa 2h, as seen in Ta- 
ble and additionally the computed Xi deviate from 
the Xi cx t-^ law (Fig. [S^a)). 

Among the other applied methods which wrongly char- 
acterize the Rl orbit as chaotic, the TDHeps scheme has 
the best performance, since Xi eventually levels off to a 
small positive value. From the results of Fig.[31[b) we see 
that the decrease of the step size r pushes the starting 
time of the level off to larger values and decreases the 
final value of Xi. So as one should expect, smaller inte- 
gration steps result in a more accurate description of the 
evolution of the orbit and deviation vectors, and leads 
to more accurate estimations of the LCEs. Nevertheless, 
the TM method is preferred over the TDHeps method 
because for the same step size r it needs less CPU time, 
and additionally it estimates more accurately the LCEs. 

For a regular orbit of the 2D Hamiltonian ((54| and a 
random choice of initial deviation vectors, the theoretical 
prediction ([TT]) for the behavior of the GALIs gives 

G2(t) cx const., G3{t)(x\, G4{t) oc \. (58) 

In Fig. m we plot the time evolution of G2, G3 and G4 
for the regular orbit Rl, when the variational equations 
are integrated by the same five numerical schemes used in 
Fig.m The results obtained by the DOP853 (Fig. g^a)), 
and the TM (Fig. IH^e)) schemes arc in accordance with 
the theoretical predictions (|58p . The GALIs computed 
by the TDHcc (Fig. IDJb)), the TDHes (Fig. g^c)) and 
the TDHeps (Fig. Hl^d)) methods follow the theoretical 
laws ([55)1 up to t « 10^ for the first two methods and 
up to < ~ 10^ for the last one. After that time the 
GALIs fall exponentially fast to zero indicating, wrongly, 
that the orbit is chaotic. This behavior is in agreement 
with the behavior of Xi obtained by these methods in 
Fig. [21 because the mLCE levels off to a positive value 
after some initial time interval, implying that the orbit 
is chaotic. The TDHeps method has again a better per- 
formance than the other two methods used to approxi- 
mate the dynamics of the TDH (|A3p . since the computed 
GALIs follow the theoretical predictions ([55)) for longer 
times, but eventually it also fails to characterize correctly 
the nature of orbit Rl. 



2. Chaotic orbits 

The computed LCEs and GALIs of the chaotic orbit 
CI are practically the same irrespectively of which of the 
previously presented methods is used for the integration 
of the variational equations. For this reason in Fig. [S] we 
present results obtained only by the DOP853 integrator. 

From the results of Fig. [5ja) we see that Xi remains 
almost constant and different from zero, having prac- 
tically the same value Xi « 4.5x 10~^ at t ~ 10® 
for all applied schemes. Thus, all used methods are 
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FIG. 4: The time evolution of G2{t) (black curves), G3{t) (gray curves) and G4(t) (light gray curves) for the regular orbit Rl 
of the Henon-Heiles system (|54|l . The variational equations are integrated by the DOP853 (a), the TDHcc (b), the TDHes (c), 
the TDHeps (d), and the TM (e) method. The plotted lines in panels (a) and (e) correspond to functions proportional to 
(dashed lines) and t^* (dotted lines). The values of r and S used in the integrations are the same as in Fig. [2] 
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FIG. 5: The time evolution of (a) Xi(t) (black curve), X2{t) (gray curve), (b) \Xi{t) + X4{t)\ (black curve), \X2{t) + Xsit)] 
(gray curve), and (c) 6*2 (f) (black curve), G3{t) (gray curve) and dit) (light gray curve) for the chaotic orbit CI of the 
Henon-Heiles system (|54p when the variational equations are integrated by the DOP853 integrator. The values of r and 5 used 
in the integrations are the same as in Fig. [2] 



able to determine correctly the chaotic nature of the or- 
bit. Since the Henon-Heiles system ((54|) is conservative, 
X2 = 0. From Fig. [D^a) we see that the finite time LCE 
X2 tends to zero, and becomes negative after t ~ 10^ 
with 1 < 10~^. At that time all the applied numeri- 
cal approaches reach their limits of applicability for the 
accurate computation of X2- The quantities \Xi -f X4I, 
IX2-I-X3I (Fig.[S]Jb)), which theoretically should be zero, 
level off after t w 10^ - 10"* to \Xi 4^ X4 1 w 4 x IQ-* and 
\X2 +^3] ~ 10^* for all used schemes. This behavior in- 
dicates that all numerical methods succeed to reveal the 
symmetric nature of the spectrum of LCEs but only up 
to four decimal digits of accuracy. Finally, the computed 
values of GALIs of orbit CI (Fig.[5{c)) show an exponen- 
tial decay to zero which is a characteristic of chaoticity. 

Fig. [5] shows the equivalence of the different numerical 
techniques in the case of the chaotic orbit CI. This is a 
clear difference with respect to the behavior of the various 
numerical schemes for the regular orbit Rl, where only 
the DOP853 and the TM methods gave similar (to each 
other) and correct results (Figs. [2] and 2]). In order to 
check if the equivalence of all methods is valid for all 



chaotic orbits we consider a weakly chaotic orbit confined 
to a thin region of the phase space at the borders of 
a small stability island (Fig. [5]). We call this orbit C2 
and its initial conditions are x ~ Q, Px ~ 0.11879, y ~ 
0.335036, py = -0.385631. 

From the results of the finite time LCEs of orbit 
C2 presented in Fig. [3 we see that both the DOP853 
(Fig. [7i;a)) and the TM (Fig. [7)^e)) methods character- 
ize orbit C2 as weakly chaotic having a small mLCE 
Xi « 4 X 10-6. The TDHcc (Fig. [7Jb)), the TDHes 
(Fig. [ZJc)) and the TDHeps (Fig. ^A)) also character- 
ize orbit C2 as chaotic but overestimate the value of 
Xi- Thus, these three methods fail to compute accu- 
rately the small value of the mLCE, with the TDHeps 
method showing once more the best performance, be- 
cause the computed value {Xi « 1.3 x lO"^) is closer to 
the real value of xi- The limitations of these three meth- 
ods are also clearly seen from the fact that the quantities 
IX1+X4I, IX2+X3I (Figs.[2i;g)-(i)) level off to larger val- 
ues with respect to the results obtained by the DOP853 
(Fig. EJf)) and the TM (Fig. E^j)) method. It is worth 
noting that the level off values of \Xi -f X4I, \X2 + Xj] 
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FIG. 6: A part of the PSS {x ^ 0, > 0) of the Henon- 
Heiles system (|54|l with H2 = 0.125, where the weakly chaotic 
orbit C2 is plotted by black dots. 

obtained for orbit C2 by the DOP853 and the TM meth- 
ods are smaller than the saturation values of the same 
quantities for the CI orbit (Fig.[5{b)). 

The results of Figs. [5] and [7] lead us to conclude that 
the DOP853 and the TM methods are able to accurately 
compute mLCEs for a larger range of xi values than the 
TDHcc, the TDHes and the TDHeps techniques. More 
specifically, our results show that the DOP853 and the 
TM schemes can evaluate xi having values at least as 
small as 10~^, while these small values definitely exceed 
the computational ability of the TDHeps method (which 
is the one with the best performance among the three 
other used methods) for the used step size r. 

The Gfc, k = 2,3,4 computed by the DOP853 
(Fig. EKa)) and the TM (Fig. l^e)) method have prac- 
tically the same behavior. Up to t w 10^, when the 
values of Xi in Figs. [7|a) and (e) start to level off de- 
viating from the Xi (x law, the GALIs follow the 
theoretical predictions ([55)1 of regular motion. Later on 
the chaotic behavior of orbit C2 becomes prominent and 
the GALIs fall exponentially to zero. The time evolution 
of GALIs computed by the TDHcc (Fig. [SKb)), the TD- 
Hes (Fig.lSIc)) and the TDHeps (Fig.[Hi;d)) method also 
indicate that the orbit is chaotic, but the exponential de- 
cay to zero starts earlier. This behavior is in accordance 
with the overestimation of orbit's chaoticity, which was 
also seen in the computation of Xi (Figs. [71[b)-(d)). 



B. Hamiltonian systems with more than two 
degrees of freedom 

Let us now apply the five different methods used in 
Sect. IVII Al to regular and chaotic orbits of the 3D and 
the 8D Hamiltonian systems ([55]) and ([55|) . In all stud- 
ied cases the computed LCEs and the GALIs have similar 
characteristics to the ones seen for the 2D system ([Ml) . 



Due to the fact that the TM, the DOP853 and the TD- 
Heps methods always exhibited the best numerical per- 
formance, we present in this section results obtained only 
by these methods for the case of regular orbits. 

In Fig. [9] we show results for the six LCEs of a regular 
orbit with initial conditions x = y = z ~ 0, px ^ 0.1, 
Py = 0.347, = (orbit R2) of the 3D system 
which was also studied in Similarly to the results 
obtained for the 2D regular orbit Rl in Fig. [2l the three 
largest finite time LCEs Xi, X2, X3 computed by the 
DOP853 (Fig. ita)) and the TM (Fig. |^c)) method, 
tend to zero following a oc i = 1,2,3, law, 
which indicates the regular nature of the orbit. These 
two methods are also able to determine the symmetric 
nature of the spectrum of LCEs, since the quantities 
\Xiit) + Xe{t)\, \X2{t) + X5{t)\ and iX^it) + X^it)] tend 
to zero (Fig.lHKd) and (f)). On the other hand, using the 
TDHeps method one would again wrongly characterize 
the orbit as chaotic because the computed Xi levels off 
at f sa 10'* to a positive value, being Xi « 1.3 x 10^'^ 
at t = 10^ (Fig. EKb)). X2 and X3 show a better con- 
vergence to zero, while the latter one becomes negative 
after t w 10^ with {X^l < 10~^. In addition, the quan- 
tity \Xi{t) + XQ{t)\ levels off to some finite value, while 
\X2{t) + X^lt)] and \X3{t) + X4{t)\ continue to approach 
zero until the end of the integration (Fig. [5Ke)). 

According to Eq. pT|) the GALIs of a regular orbit of 
the 3D Hamiltonian system ((55)1 should evolve as 

G2(t) c>c constant, G3{t) cx constant, 

(59) 

G4(t)(x^, G5(i)«i, Ge{t)(x^. 

This behavior is seen for orbit R2 in Figs. [TUT a) and 
(c) where the DOP853 and the TM method arc used 
respectively for the integration of the variational equa- 
tions. Similarly to the case of regular orbit Rl (Fig. |4]) 
the GALIs indicate that the orbit is regular. On the other 
hand, in Fig. rTOT b) where the TDHeps method is applied, 
the computed GALIs eventually show an exponential de- 
cay, wrongly suggesting that orbit R2 is chaotic. 

Finally, let us consider a particular regular orbit of 
the 8D Hamiltonian system ([55)1 which lies on a low di- 
mensional torus. In our study we impose fixed boundary 
conditions, i. e. qo{t) = q^it) = po{t) = pg{t) ~ for 
all times fix the system's parameter to /3 = 1.5, and 
consider the regular orbit with initial conditions qi = 0.1, 
Pi = 0, z = 1, 2, . . . , 8, which we call orbit R3. This orbit 
lies on a 4-dimensional torus and was also studied in ^ . 

According to the theory of GALIs developed in 10| . 
regular motion on a 4-dimensional torus implies that the 
corresponding G2, G3 and G4 remain practically con- 
stant, while the remaining indices up to Gig tend to zero 
following particular power laws (see also Fig. 4 of [l0|). 
As we can see from Fig. [Til these expected behaviors are 
weh reproduced when the DOP853 (Figs.[TTIa) and (d)) 
and the TM (Figs. [TTJc) and (f)) methods are used for 
the integration of the variational equations. On the other 
hand, the TDHeps method fails to clearly determine the 
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FIG. 7; The time evolution of Xi{t) (black curves), X2{t) (gray curves) [upper panels] and |Xi(t) + X4(t)\ (black curves), 
\X2{t) + X3{t)\ (gray curves) [lower panels] in log- log plots for the chaotic orbit C2 of the Henon-Heiles system ((54}. The 
variational equations are integrated by the DOP853 integrator ((a) and (f)), and by the TDHcc ((b) and (g)), the TDHes ((c) 
and (h)), the TDHeps ((d) and (i)) and the TM ((e) and (j)) method. The values of r and S used in the integrations are the 
same as in Fig. (2] 



DOP853 TDHcc TDHes TDHeps JM 




1234567 1234567 1234567 1234567 1234567 
loQiot log^ot log^ot loQiot loQiot 

FIG. 8: The time evolution of G2{t) (black curves), Gs,{t) (gray curves) and dit) (light gray curves) for the chaotic orbit C2 
of the Henon-Heiles system (|54|l . The variational equations are integrated by the DOP853 (a), the TDHcc (b), the TDHes (c), 
the TDHeps (d), and the TM (e) method. The values of t and 5 used in the integrations are the same as in Fig. O 



regular nature of orbit R3, as well as the dimensionality 
of the torus on which the orbit lies. From Figs. fTlT b) 
and (e) we see that the computed GALIs have a behav- 
ior similar to the one obtained by the DOP853 and the 
TM methods, which indicates the regularity of the orbit, 
but only up to t sa 10^. For t > 10^ the computed GALIs 
eventually show an exponential decay, wrongly suggest- 
ing that the orbit is chaotic. 



VIII. SUMMARY AND DISCUSSION 

We considered the problem of the accurate and fast 
integration of the variational equations of autonomous 



Hamiltonian systems. These equations govern the evo- 
lution of a deviation vector from an orbit of the system. 
The reliable determination of this evolution is necessary 
when studies of the chaotic behavior of the system are 
needed. Many chaos detection techniques, like the LCEs 
and the GALIs which we considered in our study, are 
based on the evolution of one or more deviation vectors. 

We made a detailed presentation of several numerical 
schemes for the integration of the variational equations 
and we applied them to regular and chaotic orbits of 
Hamiltonian systems with different number of degrees 
of freedom. We also investigated the efficiency of these 
methods by comparing the CPU times they need for the 
computation of the spectrum of LCEs, as well as their 
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FIG. 9: The time evolution of Xi{t), X2{t), Xz{t) (upper 
panels) and \X^{t) + X(,{t)\, |X2(t) +X5WI, \X3(t) + Xi{t)\ 
(lower panels) in log-log plots for the regular orbit R2 of the 
3D Hamiltonian system (|55|) . The variational equations are 
integrated by the DOP853 integrator ((a) and (d)), and by the 
TDHeps ((b) and (e)) and the TM ((c) and (f)) method. The 
step size is r = 0.05 for all methods. For the DOP853 method 
the parameter 5 = 10~^ is used. Dashed lines in panels (a) 
and (c) correspond to functions proportional to . 
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FIG. 10: The time evolution of Gfc(t), = 2, 3, . . . , 6 for the 
regular orbit R2 of the 3D Hamiltonian system (|55|l . The 
variational equations are integrated by the DOP853 (a), the 
TDHeps (b), and the TM (c) method. The values of r and 
5 used in the integrations are the same as in Fig. |9l The 
plotted lines in panels (a) and (c) correspond to functions 
proportional to (dashed line), (dotted line) and 
(dash-dotted line). 



ability to accurately reproduce well-known properties of 
the LCEs and the GALIs. 

The evolution of deviation vectors cannot be separated 
from the evolution of the orbit itself because the explicit 
expression of the variational equations depend on the so- 
lution of the Hamilton's equations of motion. Therefore, 
any general-purpose integration scheme for ordinary dif- 
ferential equations, like the DOP853 integrator we con- 
sidered in our study, can be used for the simultaneous 
integration of the set of equations which includes both 
the Hamilton's equations of motion and the variational 
equations. This method proved to be very reliable since 
it reproduced correctly the behavior of the LCEs and the 
GALIs for all tested orbits and systems. 



FIG. 11: The time evolution of Gk{t), k = 2,3,4,5,8 (upper 
panels) and fc = 9, 11, 13, 14, 16 (lower panels) for the regular 
orbit R3 of the 8D Hamiltonian system (|56p . The variational 
equations are integrated by the DOP853 ((a) and (d)), the 
TDHeps ((b) and (e)), and the TM ((c) and (f)) method. 
The step size is r = 0.02 for all methods. For the DOP853 
method the parameter 5 = 10~^ is used. 



When the Hamiltonian function H can be split into two 
integrable parts A and B, like H = A + B, symplectic in- 
tegrators can be used for the integration of the equations 
of motion. Symplectic integrators are known to have bet- 
ter performance than non-symplectic ones for the same 
integration time step, in terms of accuracy and required 
CPU time. In order to investigate the applicability of 
such methods for the integration of the variational equa- 
tions, we focused our study explicitly to Hamiltonians of 
the form H = A + B. In particular, we considered Hamil- 
tonians having a kinetic energy which is quadratic in the 
momenta and a potential which depends only on the posi- 
tions (Eq. (O). For such systems the two integrable parts 
A and B, are usually chosen to be the kinetic energy and 
the potential respectively. Most symplectic schemes re- 
quire the construction of symplectic maps e"^^-* (j24p and 
e"^^^ (pS)) for the solution of the integrable parts A and 
B. In our study we considered a very efheient symplectic 
method, the integrator, which has an extra degree of 
complexity with respect to most symplectic integrators, 
since it requires the explicit solution of an additional cor- 
rector term C (map e"^^*^ ([?/)) ). 

The variational equations of Hamiltonian ([5]) can be 
written as the Hamilton's equations of motion of the 
time dependent TDH (jH]), whose coefficients are defined 
by the coordinates of the orbit. Although individually 
the Hamilton's equations of motion (jB)) and the varia- 
tional equations ([7]) are equations of motion of Hamilto- 
nian functions, the system (j4ip which includes together 
both of them cannot be considered as the equations of 
motion of a new generalized Hamiltonian, and so, sym- 
plectic integrators cannot be directly used for solving it. 
In our study we applied several approaches based on sym- 
plectic techniques for the integration of the variational 
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equations. One approach we considered was the approx- 
imation of the solution of the TDH through the knowl- 
edge of the orbit's coordinates at specific times. These 
coordinates can be obtained by any symplectic or non- 
symplectic integrator, independent of the method we use 
for approximating the solution of the variational equa- 
tions. In our study we applied the integrator for this 
purpose. First we assumed the coefficients of the TDH to 
be constants for each integration step, and we integrated 
the resulting quadratic TDH by the integrator (TD- 
Hcc method) or solved it explicitely (TDHes method) 
whenever this was possible (like for example in the case 
of the Henon-Heiles system ((54)) ). An alternative way we 
also implemented was to use the integrator for inte- 
grating the time dependent TDH in an extended phase 
space (TDHeps method), using again the knowledge of 
orbit's coordinates at specific times. As an application 
of the TDHeps method we refer to the numerical study 
of the FPU problem in [slj where a leap-frog integrator 
was used for the integration of the time dependent TDH. 

The TDHcc, TDHes and TDHeps methods had a 
rather poor numerical performance as they failed in many 
cases to determine correctly the regular or chaotic na- 
ture of orbits. Our numerical results show that the com- 
puted values of the LCEs cannot become smaller than 
a small positive value, which sets a lower limit to the 
ability of these techniques to numerically determine very 
small LCEs. So, one could wrongly characterize regular 
orbits as slightly chaotic because their computed LCEs 
cannot become smaller than the above-mentioned limit, 
although their actual LCEs are zero. This happens for 
the regular orbits Rl (Fig.[2|) and R2 (Fig. [9]). Of course 
this limiting value decreases for smaller integration steps 
because the numerical schemes approximate better the 
real tangent dynamics of the system (Fig. [2Ib)). Ad- 
ditionally, one could overestimate the mLCE of chaotic 
orbits like for example in the case of the chaotic orbit 
C2 (Fig. [7]). Nevertheless these methods always required 
less CPU time than the non-symplectic DOP853 method 
for the same time step. Therefore these schemes can be 
used for some rough and fast evaluation of LCEs' charts 
but not for the detailed investigation of the dynamics or 
for the accurate computation of the LCEs and GALIs. 
We note that among these three techniques the TDHeps 
method had always the best numerical performance, al- 
though it required a bit more CPU time than the other 
two methods. 

The use of any symplectic scheme for the integration of 
the equations of motion ([5]) of the A^D Hamiltonian ([5]) 
corresponds to the repeated action of a 2A^-dimcnsional 
symplectic map S, constructed by the appropriate com- 
position of maps e^^^ ([Ml), e^^« ^ (and e^^'^ (gT]) if 
the corrector term C is used). Then, the tangent dynam- 
ics of the flow, i. e. the solution of the variational equa- 
tions (O, is described by the tangent map TS = dS/dx 
of S (some particular implementations of this approach 
for different physical problems can be found in [Sfl [13] ) ■ 

The TM method we presented in our study provides 



a simple, systematic technique to construct the tangent 
map TS for any general symplectic integration scheme 
used for the integration of the orbit, which is perfectly 
suited for practical implementations. According to this 
method, one has to substitute the 2A^-dimensional maps 
e""^-* (Ull), e^^" e^-^c' ^ needed for the symplec- 

tic integration of the equations of motion ([5]), by the 
extended 4A^-dimensional maps e^^"^^' dM]), e^^^*-' ([TT)) . 
e^'"'^^' ([53)) respectively. This procedure leads to the con- 
struction of an extended 4A^-dimensional final map com- 
posed by the 2A^-dimensional maps S and TS. In partic- 
ular, the first 2A'^ equations of this map are the equations 
of map iS", and the rest 2N equations form the tangent 
map TS. 

The TM method and the DOP853 integrator were the 
only techniques that succeeded in computing correctly 
the LCEs and the GALIs for all studied cases. Among 
them, the TM method required less CPU time for the 
same integration step size. Another advantage of the 
TM method over the DOP853 integrator is that its ap- 
plication with larger time steps reduces the needed CPU 
time, keeps the accuracy to acceptable levels, and pro- 
duce more reliable results than the DOP853 integrator. 

In conclusion, the TM method proved to be the most 
efhcient one among all tested methods, since it required 
the least CPU time for the computation of the spectrum 
of LCEs and reproduced very accurately the behavior of 
the LCEs and GALIs. Therefore, whenever the studied 
Hamiltonian can be split into two integrable parts, so 
that it can be integrated by symplectic integrators, the 
TM method should be preferred over other symplectic or 
non-symplectic integration schemes. 

Although we considered in our study applications of 
the TM method to Hamiltonian systems of relatively low 
dimensionality (systems having up to eight degrees of 
freedom), the method is expected to be also very efficient 
for higher-dimensional systems. Symplectic integrators 
have already been applied successfully for the accurate 
integration of motion in multi-dimensional systems which 
are related for example, to problems of astronomical in- 
terest (e. g. [131), of molecular dynamics (e. g. (29l. [ssj) 
and dynamics of nonlinear lattices (e. g. [23] ) . Using the 
TM method these symplectic integration schemes can be 
extended to integrate also the corresponding variational 
equations. This is a problem of great practical impor- 
tance, which we plan to address in a future publication. 

As a final remark, we note that all the presented meth- 
ods require the knowledge of the analytic expression of 
matrix Dy((f(t)) ^ (or of matrix A(t) @ in the case 
of a general dynamical system). If the variational equa- 
tions cannot be written explicitly, possibly due to the 
complicated form of the studied dynamical system, the 
analytical derivation of these matrices is not possible and 
their elements could be estimated numerically, introduc- 
ing an additional error to the solution of the variational 
equations. An approach that could be followed in such 
cases is the approximation of the solution of the varia- 
tional equations by the difference of two orbits initially 
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located very close to each other (see [39| for some partic- 
ular applications of this approach). This is the so-called 
two-particle method, which was introduced in [l^ and is 
mainly used for the evaluation of the mLCE. It was re- 
alized almost immediately after the introduction of this 
technique that this approach is less efficient and reliable 
than the actual integration of the variational equations 
[33} (whenever, of course, this integration is possible). 
For this reason we did not include this approach in our 
study. 
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1. Symplectic integration of the equations of 
motion 

The Hcnon-Heiles Hamiltonian can be split into 
two parts H2 — A + B, according to equation ((23)) . with 

A = lipl+pl), 

(A4) 

D 1 / 2 I 2\ I 2 1 3 

B = -{x +y )+x y - -y . 

As it was explained in section |Vl this separation is con- 
venient for the application of symplectic schemes for the 
integration of equations (jAip . since Hamiltonians A and 
B are integrable. The maps e'^^'^ ([24]), e'^-^^ ((251), which 
propagate the set of initial conditions (x, y-,Px,Py) a-t time 
t, to their final values {x' , y' ,p'x-,p'y) at time t-\- t are 



X' 


= X+PxT 


y' 


= y+PyT 


p'. 


Px 


P'y 


Py 



(A5) 



Appendix A: Analytical expressions for the 
integration of the Henon-Heiles system 

We present here the explicit expressions of the vari- 
ous integration schemes for the 2D Hcnon-Heiles system, 
whose Hamiltonian function (|54p is of the form ([5]) with 
q = {x,y), p = {px,Py). The Hamilton's equations of 
motion ^ are 



X 

y 



= Px - x{l + 2y)T 
= Py + (y^ - - y)T 

The corrector term ([26]) is 

C = {B, {B, A}} - [x + 2xyf + {x^ - y^ + y) 
and the corresponding map e"^^^ (P?]) 



(A6) 



(A7) 



X 

y 

Px 

Py 



Px 
Py 



2xy 



The variational equations ^ of the system are 

6x — 6px 

6y = 6py 

5p^ = -~{l + 2y)Sx~2xSy ' 

6py = -2xSx + {-1 + 2y)Sy 



^tLc . ) y = V 

(Al) ■ \ p', ^ Px-2x{l + 2x'' + &y + 2y'')T 

p'y = Py - 2{y - 3y2 + 2y^ + 3x^ + 2x^y)T 

(AS) 



2. Integration of the variational equations 

We derive now for the particular case of the Henon- 
(A2) Hciles system the analytical expressions of the various 
numerical schemes presented in section IVII for the inte- 
gration of the variational equations. 



while the corresponding TDH ^ takes the form 



Hvh{Sx, Sy, Spx, Spy] = 2 i^Pl + ^pD + 
+ i { [1 + 2yit)] Sx^ + [1 - 2y{t)] Sy^ + 2 [2x{t)] SxSy} 



(A3) 



a. Diagonal form of the TDH lfA3\) with constant 
coefficients 

Inserting the values x(ti) = Xi, y{ti) = yt at a specific 
time ti in the functional form of the TDH (|A3|) . Hyu 
becomes a quadratic 2D Hamiltonian with constant co- 
efficients. The equations of motion of this Hamiltonian 
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are solved immediately if Xi 
formation 



0. For Xi ^ Q the trans- with 



5x 
5y 



SX 
SY 



5px 

5py 



SPx 
6Py 



(A9) 



(AlO) 



gives HyniSx, Sy, 5px, 6py\ ti) the diagonal form 
Hvhd{SX, 6Y, SPx,SPy) = i {SP^ + 6P^) 



l + 2Jx^^+yf SX^+[l-2Jx^^+yf SY 



(All) 

The columns of matrix T are the eigenvectors of matrix 



T>l{q{U)) = -Dlix,,y,) 



1 + 2yi 2x, 
2xi 1 - 2yi 



, (A12) 



and Ai,2 = l±2y^xf + yf the corresponding eigenvalues. 



6. Symplectic integration of the TDH in an extended 
phase space 



Considering the TDH (jASp time dependent Hamil- 
tonian, we can transform it to a time independent one 
having time t as an additional generalized position by the 
procedure presented in section IVI B 21 The 3D Hamilto- 
nian takes the form 

HvH{Sx,6y,t,6px,Spy,pt) = ^ {6pl + Spl) +pt 



+ - { [1 + 2y{t)] Sx" + [1 - 2y{t)] Sy'' + 2 [2x{t)] SxSy) , 

(A13) 

with pt being the conjugate momentum of coordinate t. 
HvH can be split into two integrable parts ([33]) 

A{Spx,Spy,pt) = ^(Spl + Spl) +pt, 



B{Sx,Sy,t) - -{[l + 2y{t)]5x' + [l-2yit)]Sy'+ 



2 [2x{t)] SxSy} , 



(A14) 



so that its equations of motion can be integrated by any 
symplectic integration method in order to obtain the time 
evolution of variations Sx, Sy, Spx, Spy. The maps e"^^-? 
(pi]) . e'^^B (|35|) (neglecting the equations for pt) are 



( Sx' 


= Sx + SpxT 


W 


= Sy + SpyT 




= t + T 


Sp'x 


= Spx 


I Sp'y 


^ SPy 



(A15) 



= Sx 
= Sy 
= t 

= Spx - {[1 + 2y{t)] Sx + 2x{t)Sy} t 
= Spy + {-2x{t)Sx + [~1 + 2y{t)] Sy] r 

(A16) 

The corrector term C (l36l) is 




C^[Sx + 2x{t)5y + 2y{t)Sxf+[Sy + 2x{t)Sx - 2y{t)Syf , 



and the corresponding map e'' c- (|37| 



(A17) 



Sx' = Sx 
Sy' = Sy 
t' = t 

Spx - 2 {'ix{t)Sy+ 
+ [4x'(t) + {l + 2y{t)f 
Spy -~2{4x{t)Sx+ 
+ [4x'{t) + {l-2yit)f] Sy} 



Sp'x 
Sp'y 



Sx\ 



(A18) 



The tangent map method 



According to the TM method presented in section lVI CI 
equations (jAip and (IA2[) form a set of equations which 
defines the act of the differential operator Lhv on vector 
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u = {x,y,px,Py,Sx,Sy,Spx,5py) (equations (|4T|) '). This 
set of equations is split into two integrable sets 



y 

Px 
Py 

Sx 
Sy 
SPx 
^Pv 



Px 
Py 




Spx 

Spy 






du 
IE 



Lavu, 



(A19) 



y 

P'x 



y 

Px - x{l + 2y)T 
Py + (y^ - a;^ - y)T 
5x 
Sy 

Spx - [(1 + 2y)6x + 2x6y] r 
Spy + [-2xdx + (-1 + 2y)5y] r 

(A22) 

while the map e^^^^' ((53)) of the corrector function C 
(Ell) is 



Py 

Sx' 

Sy' 

Sp'x 

ysp'y 



y 

Px 

Py 
Sx 
Sy 

SPx 

SPy 



2xy 

,^2 



y 



y 




-{l + 2y)Sx~ 
-2xSx + (-1 



2xSy 
h2y)Sy J 



du 



(A20) 

which define the act of operators Lav (|44)) and Lbv (|45)) 
respectively. Then, maps e^'^-^'^ dM]) and e'^^^^ (07]) are 



r x' 


= x + 


PxT 


2/' 


= y + 


PyT 


px' 


= Px 




py' 


= Py 




Sx' 


= Sx - 


f SpxT 


Sy' 


= Sy- 


1- SpyT 


5p'x 


- Spx 




[Sp'y 


= Spy 





p'x 

Py 
Sx' 

Sy' 

6p'x 

Sp', 



2x{l + 2x^ 
2(y-3y2 + 



f 6y 
2y'- 



(A21) 



x 

y 

Px 

Py 
Sx 

Sy 

Spx - 2 [(1 + 6x2 ^ 2y^ + 
+2x{3 + 2y)Sy] r 
Spy - 2 [2x{i + 2y)Sx+ 
+ (1 + 2a;2 +6y2 - Qy)Sy] 



2a;2y)r 



(!>y)Sx- 



(A23) 
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